home *** CD-ROM | disk | FTP | other *** search
Text File | 1991-10-05 | 4.7 KB | 164 lines | [TEXT/MPS ] |
- UNIT RandomNumbers;
-
- { This unit implements a "minimal standard" random number
- generator using 32-bit integer arithmetic. This is not
- the "best" random number generator possible, but it is
- simple, quick, and produces numbers which are "random"
- enough for most purposes.
-
- It is based on the "parametric multiplicative linear con-
- gruential algorithm" which was originally proposed by D.
- H. Lehmer in 1951. This algorithm generates a sequence
- of integers z[1], z[2], z[3]... by repeatedly carrying
- out the following calculation:
-
- z[n+1] = (a * z[n]) mod m
-
- where "a" and "m" are suitable constants. The successive
- "z"s are stored in the variable "seed" in this unit.
- They all lie in the range [0..m-1], so dividing them by
- "m" gives real numbers which lie in the range [0.0..1.0].
-
- In 1969, Lewis, Goodman and Miller proposed the choice of
- a = 16807 and m = 2**31 - 1 = 2147483647 as particularly
- suited for a machine with a 32-bit word length (in their
- case, the IBM System/360). The main problem is that the
- product a * z[n] can be more than 32 bits long, leading
- to integer overflow on many systems (including the Mac).
-
- In 1979, G. L. Schrage devised an implementation of the
- LG&M method which is guaranteed not to produce integer
- overflow on any system with a 32-bit word length. It is
- this implementation which is used here, in procedure
- "UpdateSeed".
-
- Source: Stephen K. Park and Keith W. Miller, "Random
- Number Generators: Good Ones Are Hard to Find", Communi-
- cations of the ACM, vol. 31, p. 1192 (October 1988).
-
- This unit was written in MPW Pascal, v3.2.
-
- Jon Bell
- Dept. of Physics & Computer Science
- Presbyterian College
- Clinton SC 29325
- CompuServe: #70441,353
- October 1991 }
-
- INTERFACE
-
- {----------------------------------------------------------}
-
- procedure InitRandomSeed (newSeed : longint);
-
- { Initializes the random number seed to "newSeed". You
- must call this procedure once, at the beginning of your
- program, before you use any of the following functions.
- As far as randomness is concerned, it makes no difference
- what value you use for "newSeed", so long as it isn't
- zero. Using different seeds merely gives you different
- sequences of "random" numbers. Using the same seed each
- time you run the program gives you the same sequence of
- "random" numbers each time, which may be useful for
- debugging. }
-
- {----------------------------------------------------------}
-
- function RandomSeed : longint;
-
- { Returns the current value of the random number seed. }
-
- {----------------------------------------------------------}
-
- function RandomReal : extended;
-
- { Returns a real number, uniformly "randomly" selected
- from the interval (0.0,1.0). Note that this is an
- "open" interval, that is, it does _not_ include 0.0 or
- 1.0. }
-
- {----------------------------------------------------------}
-
- function RandomInteger (max : longint) : longint;
-
- { Returns an integer, uniformly "randomly" selected from
- the interval [1,max]. Note that this is a "closed"
- interval, that is, it _does_ include 1 and max. }
-
- {----------------------------------------------------------}
- {----------------------------------------------------------}
-
- IMPLEMENTATION
-
- const
- a = 16807;
- m = 2147483647;
- q = 127773; { m div a }
- r = 2836; { m mod a }
-
- { NOTE: Other possible values for these constants, which
- may give "better" results, are:
-
- a = 48271, q = 44488, r = 3399 or
- a = 69621, q = 30845, r = 23902,
-
- keeping m = 2147483647. }
-
- var
- seed : longint;
-
- {----------------------------------------------------------}
-
- procedure InitRandomSeed (newSeed : longint);
-
- begin
- seed := newSeed
- end;
-
- {----------------------------------------------------------}
-
- function RandomSeed : longint;
-
- begin
- RandomSeed := seed
- end;
-
- {----------------------------------------------------------}
-
- {private} procedure UpdateSeed;
-
- var
- lo, hi, test : longint;
-
- begin
- hi := seed div q;
- lo := seed mod q;
- test := a * lo - r * hi;
- if test > 0 then
- seed := test
- else
- seed := test + m
- end;
-
- {----------------------------------------------------------}
-
- function RandomReal : extended;
-
- begin
- UpdateSeed;
- RandomReal := seed / m
- end;
-
- {----------------------------------------------------------}
-
- function RandomInteger (max : longint) : longint;
-
- begin
- UpdateSeed;
- RandomInteger := (seed mod max) + 1
- end;
-
- {----------------------------------------------------------}
-
- END.
-